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A HYSTERESIS MODEL FOR PIEZOCERAMIC MATERIALS 


RALPH C. SMITH* AND ZOUBEIDA OUNAIES+ 

Abstract. This paper addresses the modeling of nonlinear constitutive relations and hysteresis inherent 
to piezoceramic materials at moderate to high drive levels. Such models are necessary to realize the full 
potential of the materials in high performance control applications, and a necessary prerequisite is the 
development of techniques which permit control implementation. The approach employed here is based on 
the quantification of reversible and irreversible domain wall motion in response to applied electric fields. A 
comparison with experimental data illustrates that because the resulting ODE model is physic-based, it can 
be employed for both characterization and prediction of polarization levels throughout the range of actuator 
operation. Finally, the ODE formulation is amenable to inversion which facilitates the development of an 
inverse compensator for linear control design. 
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1. Introduction. Piezoceramic materials are being employed in an increasing number of structural 
and structural acoustic systems due to the fact that they are lightweight, compact, relatively inexpensive 
and provide the capability for both sensing and actuating. These advantages are bolstered by the property 
that they exhibit nearly linear dynamics and minimal hysteresis at low drive levels. The restriction of 
the materials to the quasi-linear range, however, limits their output to levels which are a fraction of those 
attainable at moderate to high drive levels. The tradeoff for improved performance is the manifestation of 
the hysteresis inherent to the materials as illustrated in Figure 1 with data from a PZT 5 A wafer. To utilize 
the full capability of piezoceramic actuators, it is necessary to model the nonlinear constitutive relations and 
hysteresis in a manner which is conducive to both system design and control. In this paper, we consider a 
model based on the domain properties of the materials. 

The model is based on the theory presented in [9, 10] where it was illustrated for the relaxor ferroelectric 
material PMN-PT-BT employed in its ferroelectric regime. In light of this theory, the modeling of hysteresis 
in piezoceramic materials is considered in two steps. In the first, the hysteresis-free, or anhysteretic, relation 
between the applied field and resulting polarization is quantified through a classical application of Boltzmann 
statistics. This provides a constitutive relation which is applicable at low drive levels when hysteresis is 
negligible. Hysteresis is then modeled through the characterization of reversible and irreversible motion of 
domain wall pinned at inclusions in the material. The combination of the components provides an ODE model 
which incorporates the nonlinear constitutive relations and hysteresis observed in various PZT compounds 
at high drive levels. 
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To place this model in the context of existing hysteresis models for piezoceramic materials, we note 
that models can be roughly categorized as microscopic, macroscopic and semi-macroscopic in nature. The 
microscopic theories focus on the quantification of energy mechanisms at the lattice or domain level. As 
illustrated in [7], such approaches can provide a high degree of detail but are limited to simple stoichiometries 
and are currently difficult to implement in control design due to the large number of required parameters. 
Macroscopic theories lie at the opposite end of the spectrum and are advantageous when the underlying 
physical mechanisms are poorly understood or difficult to implement. This category includes models based 
on anhysteretic curves shifted to coincide with remanence polarizations [13] as well as the Preisach models 
which have been developed to quantify hysteresis in piezoceramic materials [1, 2, 3]. While advantageous 
in control design, these models often require a large number of nonphysical parameters which can limit 
their utility in systems which exhibit transient dynamics or slowly varying operating conditions. The model 
employed here is semi-macroscopic in the sense that it is based on energy" relations for the material but 
employs macroscopic averages to obtain a small number (five) of effective parameters for certain material 
properties. Due to the physical nature of the parameters, they are easily obtained when specifying the 
model for particular applications and can be efficiently updated to accommodate slowly changing operating 
conditions. 

The model is summarized in Section 2 arid validated through a comparison with data from a PZT 5 A 
wafer in Section 3. The experimental results illustrate both the accuracy of the model and its capability for 
predicting the polarization due to changing drive levels. 



-2 0 2 
Electric Field (MV/m) 


Figure 1 . Hysteresis observed in a PZT 5A wafer in response to a 1600 V input. 


2. Hysteresis Model. We summarize here the theory presented in [9, 10] for quantifying hysteresis in 
ferroelectric materials. This development is considered in two steps. In the first, the anhysteretic polarization 
is modeled through classical Boltzmann statistics. This provides the minimal energy state which would be 
attained by the material in the absence of pinning sites which restrict domain wall movement. Hysteresis 
is incorporated in the second step through the quantification of the energy required to bend and translate 
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domain walls. This respectively provides reversible and irreversible components to the polarization which 
can be summed to obtain the total polarization for the material. 

2.1. Anhysteretic Polarization. The anhysteretic polarization curve is that which would be obtained 
in an idealized material which is devoid of inclusions or imperfections. As illustrated in Figure 2, the 
anhysteretic has a burst region near the origin, where dipole rotation toward preferred ionic configurations 
produces large changes in polarization, and exhibits saturation at high field levels when charge distributions 
prohibit further changes. In contrast to the measured polarization curves at high drive levels, the anhysteretic 
is single- valued and reversible. 

To quantify the dependence of the anhysteretic polarization on field and coupling effects, w r e consider the 
effective field acting on moments in the material. Under the assumption of fixed temperatures, no applied 
stresses and quasi-static operating conditions, the effective field is taken to be 


( 2 . 1 ) 


E e = E + aP 


where E and P respectively denote the applied field and resulting polarization, and the term aP quantifies 
the field due to interdomain coupling. The parameter a = Eq/P s , where Eq and P s respectively denote 
a scaling electric field and the saturation polarization, is typically estimated through a least squares fit 
to data. We note that throughout this development, P s denotes the theoretical saturation value beyond 
which polar interactions prevent further increases in polarization. This definition is in accordance with the 
corresponding definition for the saturation magnetization (e.g., see [5]) and should not be confused with the 
definition employed in several texts (e.g., [4, page 36]) for the saturation value obtained by extending the 
slope at tip reversal back through the vertical axis. 

As detailed in [9, 10], Boltzmann statistics is then employed to quantify the probability of dipoles 
occupying certain energy states. The anhysteretic models differ according to the assumptions regarding the 
possible dipole orientations and dielectric behavior of the material. Under the assumption that the material 
is isotropic and the orientation of cells can be in any direction, the balance of the thermal and electrostatic 
energies through a Boltzmann probability relation yields the Langevin equation 


( 2 . 2 ) 


P, 



a 

~E e 


for the anhysteretic polarization. The constant a is specified by a = where T and T c are the operating 
and Curie temperatures for the material. Because the scaling field Eq is unknown, this constant must also 
be estimated for a given material. 

A second model which is commonly employed in ferroelectric applications is the Ising spin relation 
(2.3) P an = P s tanh 


which results from the assumption that dipole moments can only be oriented in the direction of the field 
or opposite to it. As noted in Figure 2 where the two models are compared, the restrictions on possible 
moment orientation causes the Ising spin relation to saturate more quickly than the Langevin model. Further 
discussion concerning the merits of the two relations as well as details regarding their derivation are provided 
in [9, 10]. 
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Figure 2. Ising spin and Langevin models for the anhysteretic polarization. 


2.2. Domain Wall Model. The relations (2.2) or (2.3) can be employed at low drive levels but are 
inappropriate at moderate to high drive levels since they do not incorporate the hysteresis inherent to the 
materials. As detailed in [10] and included references, sigmoidal hysteresis of the type depicted in Figure 1 
is typically attributed to the energy required to translate domain walls across pinning sites in the material. 
The pinning of domain walls by defects or inclusions in the material is due to the reduction in energy 
which occurs at these sites. The equilibrium position for the domain wall occurs when the potential energy 
is minimized. At low field levels, the walls remain close to the equilibrium position and the motion is 
reversible. From an energy perspective, the variations are not sufficient to cross a barrier in the potential 
well. The motion becomes irreversible when sufficient energy is provided to cross the potential barrier. 
Physically, this can occur when the domain wall intersects a remote pinning site and is the mechanism 
underlying domain wall translations. The resulting irreversible polarization Pj rr and reversible polarization 
P rev are then summed to obtain the total polarization. This approach follows that employed by Jiles and 
Atherton in their corresponding hysteresis model for ferromagnetic materials [6]. 

To quantify the irreversible polarization, it is noted in [9, 10] that the polarization level for a given 
effective field can be expressed as that for the ideal case minus losses required to break pinning sites. This 
yields the relation 

(2.4) • Pirr = Pan ~ Sk ^ 

where the parameter <5 = sign(dE) ensures that the energy required to break pinning sites always opposes 
changes in polarization. The parameter k is defined by k = -4^ where n denotes the average density of 
pinning sites, (£*) is the average energy for 180° walls and p is an average dipole moment. Because the 
density and energy of individual pinning sites are unknown, the parameter k must be estimated for a given 
material. 

The formulation of (2.4) in terms of the applied field E yields the differential equation 


dPirr 

dE 


Pan ~ Pit 


Sk Ot(Pan Pirr) 
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specifying the irreversible polarization. As discussed in [9, 10], while this expression is adequate in most 
regimes, it can yield nonphysical solutions when the field is reversed from saturation for materials which 
exhibit significant hysteresis and are driven at high levels. The enforcement of solely reversible polarization 
changes in this regime eliminates this discrepancy and yields the relation 

fry r ^ dPfrr T P 3 71 P ITT 

{ dE ~ kS - a (P on - P irr ) 

where 

~ f 1 , {dE > 0 and P > P an } or {dE < 0 and P < P an } 

[ 0 , otherwise. 

The second component of the polarization is the' reversible polarization which models the effects of 
domain wall bending. To first approximation, this is mocieled by the relation 

(2.6) P rev = c(P an — Pirr) 

where c is a parameter which must be estimated for the specific application (see [9, 10]). 

The total polarization is then given by 

P = Prev T Pirr 


or equivalently, 

(2.7) P = cP an + (1 - c)P irr . 

To implement the model, the effective field for a given field and irreversible polarization level is computed 
using (2.1). This effective field value is then employed in either (2.2) or (2.3) to compute the corresponding 
anhysteretic polarization. The subsequent irreversible polarization is determined by numerically integrating 
(2.5). The total polarization is then specified by (2.7). Various aspects of the model are illustrated in the 
next section through a comparison with experimental data. 

3. Model Validation. To illustrate the performance of the model, we consider the characterization of 
the polarization generated by a cylindrical PZT 5A wafer. The wafer had a diameter of 1 in and was 10 
mils thick. To maintain quasistatic operating conditions, the frequency of the input electric field was taken 
to be 200 mHz. Because the material was initially polarized prior to the hysteresis measurements, the first 
cycle was omitted because it exhibited transient behavior. Three complete steady state cycles were then 
measured for input voltages ranging from 600 V to 1600 V. The corresponding field inputs to the model were 
determined using the relation E = V/d where d— 10 mils is the thickness of the wafer. 

Two strategies were considered to illustrate the capability of the model to quantify the inherently hys- 
teretic relation between the input field and generated polarization at various drive levels. In the first, the 
parameters a,a,c, k and P s were obtained through a least squares fit to polarization data collected with a 
peak input voltage of 600 V. The model, with these parameters, was then used to predict the polarization 
at drive levels up through 1600 V. The second option differed in that parameters were obtained by fitting 
the 1600 V data and the model was then used to predict the polarization behavior at lower drive levels. 

The model behavior and corresponding data for the first strategy are illustrated in Figure 3. The 
parameter values a = 3.7 x 10 6 Vm/C, a = 4.1 xlO 5 C/m 2 , c = 0.35, k = 1.8 xlO 6 C/m 2 and P s = 0.49 C/m 2 
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were obtained through a least squares fit to the 600 V data with the Langevin relation used to model the 
anhysteretic polarization. The model predictions at 800 V, 1000 V and 1600 V were then compared with 
the corresponding data at these drive levels. Throughout the range of drive levels, the model accurately 
predicts the measured polarization. We note that the capability of the model for predicting the polarization 
at various drive levels is due to the fact that it is based on energy principles. 

The identification of parameters at 1600 V produced model fits that were slightly more accurate at that 
level but were less accurate when predicting at 600 V. For this data set, the 600 V data appears to provide 
more information relevant to identifying the parameters and hence gives the model slightly better predictive 
capabilities than identification at high drive levels. For general PZT applications, however, determination 
of parameters and comparison of model performance at both low and high drive levels will indicate which 
parameters provide the best overall fit. An advantage of this approach is that once the parameters are 
determined, the model can be used to quantify the polarization for a large range of drive levels simply by 
specifying the input voltage or field. 

4. Concluding Remarks. This paper summarizes a hysteresis model for piczoceramic actuators which 
is based on domain wall theory for ferroelectric materials [9, 10]. The theory quantifies the hysteresis inherent 
to the materials through the consideration of the energy required to bend and translate domain walls which 
are pinned at defects or inclusions in the material. This provides reversible and irreversible polarization 
components whose sum represents the polarization due to an applied field. Characterization in this manner 
provides the model with the capability for specifying the polarization at a variety of input field levels with 
one set of model parameters. The flexibility of the model is further augmented by the small number (five) 
of required parameters and the physical nature of the parameters (e.g., P s is often known a priori or can be 
directly obtained from the data). 

The performance of the model was illustrated through a comparison with data from a PZT 5A wafer. 
The resulting model behavior was accurate both at the drive level used to estimate the parameters as well 
as input fields up to more than double that magnitude. This illustrated the flexibility of the model and its 
capability for prediction throughout the drive range of the actuator. 

Finally, the ODE nature of the model makes it amenable to inversion through the consideration of a 
complementary ODE in a manner analogous to that described in [8]. This facilitates the construction of an 
inverse compensator which can be used for linear control design [11, 12]. The application of these techniques 
for linear control implementation for piezoceramic actuators which exhibit hysteresis is under investigation. 


Acknowledgements. The research of R.C.S. was supported in part by the Air Force Office of Scientific 
Research under the grant AFOSR-F49620-98-1-0180. 
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Figure 3. Model fit to 600 V data and model predictions at 800 V, 1000 V and 1600 V with the parameter 
choices a = 3.7 x 10 6 Vm/C, a = 4.1 x 10 5 C/m 2 , c = 0.35, k = 1.8 x 10 6 C/m 2 and P s = 0.49 C/m 2 ; Model 
( — — ), Data ( ). 
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